! This represents  n_modulation + N_spin = integer resonance
!
! This is a program which uses a version of the ALS represented by a flat file
! The sextupoles are turned off
! The symmetry is broken naturally by modulating a single element, namely BEND1
! The entire acceleration is faked a la Lucio (God help us)
!

program ptc_geometry
use madx_ptc_module
use pointer_lattice
use c_TPSA
implicit none

character*48 :: command_gino
type(probe) xs0,xs1,XST
type(probe_8) xs
type(layout), pointer :: als
real(dp) closed(6),BETA_XX,BETA_XY,theta0,N0(3),DN0_DX(3),DN2_DX2(3),norm,dx,X0(6),XN(LNV)
type(real_8) ray8(6),theta
type(damap) id,one_turn,a
type(normalform) nf
INTEGER MF,I,N,k,mft,pos
TYPE(FIBRE),POINTER:: P
TYPE(integration_node),POINTER:: t,t1,t2
type(internal_state) state
type(damapspin) id_s,one_turn_S,a_S,A_f,A_spin,A_l,A_nl,NORMALSPIN,M_DS,R_DS
type(normal_spin) nf_S
TYPE(spinor_8) N_AXIS
TYPE(spinor) N_AXIS_ave
TYPE(res_spinor_8) N_AXIS_res
real(dp) dr,r,ag,ag0,gamma0,AR,AI,bR,bI,cR,cI,at,dre,di,art,aitt,fac
integer spinres,number_of_ac_plane
real(dp) x(6),circ,xt(lnv),prec 
complex(dp) xc(6) , bet, dtune
logical(lp) mis
TYPE(work) E_0
REAL(DP) a13,a11,rate
type(c_ray) cray,xcn
!-----------------------------------

type(c_damap) c_map,c_spin0
type(c_normal_form) c_n
TYPE(c_spinor) C_N0
    interface
       subroutine build_ALS(ALS,MIS)
         use madx_ptc_module
         use pointer_lattice
         implicit none
         type(layout), target :: ALS
         logical(lp) mis
       end subroutine build_ALS
    end interface



Lmax = 10.d0
use_info = .true.
spinres=1
!write(6,*) " spin resonance 1 or -1 "
!read(5,*) spinres
!== user stuff : one layout necessary before starting GUI
!call ptc_ini
call ptc_ini_no_append
ALS=>m_u%start

!Write(6,*) " small misalignments and rotations in lattice ? input t or f "
!read(5,*) mis

!call build_ALS(ALS,mis) 

!call ptc_ini_no_append

CALL  READ_AND_APPEND_VIRGIN_general(M_U,"ALS_FLAT_fast.DAT")
m_U%start%name="ALS "
write(6,*) "Making an node layout"
CALL MAKE_NODE_LAYOUT(m_U%start)

als=>m_U%start
p=>als%start


!!!! circ is the circumference of the ring !!!! 
call get_length(als,circ)

!!!! modulate.txt sets the magnet BEND1 as a modulated magnet !!!! 
call read_ptc_command77("modulate.txt")

!


 !  
ag0=  1.158630183290000D-003 
p=>als%start
E_0=p
write(6,*) 1.d0/e_0%gamma0I,p%ag/e_0%gamma0I
gamma0=1.d0/e_0%gamma0I
 

!!!!!  Make the lattice linear by zeroing all sextupoles !!!!!! 
!!!!!  Set all magnets spin to g=ag0 for fake electron tracking!!!!!! 
p=>als%start
do i=1,als%n
 p%ag = ag0  
 IF(P%MAG%P%NMUL>=3) THEN
  CALL ADD(P,3,0,0.D0)
 ENDIF
 p=>p%next
enddo
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

closed=0.d0
state=default0+nocavity0  

    CALL FIND_ORBIT(ALS,CLOSED,1,STATE,c_1d_5)
 
call kanalnummer(mf,"spin_result_barber_rf.txt")

state=default0+only_4d0+SPIN0
XS0=CLOSED
    CALL FIND_ORBIT_probe_spin(ALS,XS0,STATE,c_1d_5,FIBRE1=1,theta0=theta0)

WRITE(6,*) XS0%X
WRITE(6,*) XS0%S(1)%X
WRITE(6,*) THETA0/twopi
N0=XS0%S(1)%X
write(6,*) " n0 "
write(6,*) n0

state=default0+only_4d0+SPIN0+modulation0


!!!!  The last dimension of phase space are a pseudo-clock modulation the bend
!!!!! (Pseudo because it uses path length along some ideal orbit to advance
!!!! it is carried along and when the bend is traverse, the field is modulated
!!!! according to the x-component of the clock 
!!!! set a modulation clock !!!!!!
xs0%ac%om=0.401073997029287d0*twopi/circ
xs0%ac%x(1)=1.d0
xs0%ac%x(2)=0.d0
write(6,*) " Modulation tune =",circ*xs0%ac%om/twopi

!!!! xs0%ac%x(1) and xs0%ac%x(2) will be the (q,p) of the modulating clock
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
call init_all(STATE,3,0)

call alloc(c_map)
call alloc(c_n)
call alloc(c_spin0)
CALL ALLOC(C_N0)
!print_spin=.true.

    call alloc(xs)
    call alloc(ray8)
    call alloc(theta)
    call alloc(N_AXIS_res)
    call alloc(N_AXIS)
    call alloc(id_S,one_turn_S,A_S,A_f,A_spin,A_l,A_nl,NORMALSPIN)
    call alloc(nf_S)

!!!! Polymorphic probe is created in the usual manner   
   ID_S=1
   
   XS=XS0+ID_S     !
   
!!!! get spin polymorphic probe after one turn   
CALL TRACK_PROBE(ALS,XS,STATE,FIBRE1=1)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 

!!!! Copy probe_8 into a complex damap because we use the complex package for normal forms
c_map=XS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
 
!!!! one resoanance is left in the map 
c_n%nres=1
c_n%m(3,1)=1   !  the third plane here is modulation see above
c_n%ms(1)=1    ! spin tune
!    so the resonance is   n_modulation + N_spin = integer
call c_normal(c_map,c_n,dospin=my_true)
prec=1.d-10
! printing all the tunes
write(6,*) c_n%tune(1:c_%nd)

write(6,*) " spin ",c_n%spin_tune

!  Here I put A_t in phasor's basis!!!    
c_n%A_t=c_n%A_t*from_phasor()   

! m = As a n a^-1 As^-1
c_map=c_n%A_t**(-1)*c_n%As**(-1)*c_map*c_n%As*c_n%A_t
 
write(mf,*) " Normalised map "
call  print(c_map,mf,DEPS=prec)

!!! C_spin0 is a complex map
c_spin0=1   !!! identity

!!!! Only contains linear part of spin, 0-zeroth order matrix
C_SPIN0%S=c_map%S.SUB.0

write(mf,*) " Linear part of the spin only  "
call  print(C_SPIN0,mf,DEPS=prec)

! Here the map is factorised as map :   exp(theta0 L_y) o ( r, exp(a _y + resonance))
! C_SPIN0 =( r, exp(a _y + resonance))

C_SPIN0=C_SPIN0**(-1)*C_MAP

write(mf,*) " "
write(mf,*) " Here the map is factorised as map :   "
write(mf,*) "  exp(theta0 L_y) o ( r, exp(a _y + resonance)) "
write(mf,*) "   ( r, exp(a _y + resonance)) is printed"
write(mf,*) "   "

call  print(C_SPIN0,mf,DEPS=prec)


! a _y + resonance = C_n0
write(mf,*) " "
write(mf,*) "   a _y + resonance = C_n0  "
write(mf,*) " "
          call c_find_n0_da(c_spin0%s,C_N0)
          call c_n0_to_nr(C_N0,C_N0)
write(mf,*) " C_n0 is  "

Write(mf,*) " C_n0%v(1)  < =   C_n0%v(1)-i_*C_n0%v(3) ! coefficient of  1/2(L_x + i L_z) "
Write(mf,*) " C_n0%v(2)  < =   C_n0%v(2)"
Write(mf,*) " C_n0%v(3)  < =   C_n0%v(1)+i_*C_n0%v(3) ! coefficient of  1/2(L_x - i L_z) "
write(mf,*) " "
CALL PRINT(C_N0,mf,PREC)

!!! Here starts some evaluation of the various operators

write(mf,*) " "
write(mf,*) " Here I set the initial ray of the tracking"
write(mf,*) " "

fac=1.7d0
Xc=0.D0
Xc(1)=0.14D0*fac
Xc(3)=0.01D0*fac
Xc(5)=1.d0
XcN%x=0.D0

write(6,*)
cray=xc
write(mf,'(6(1x,g12.5))') xc(1:6)

write(mf,*) " "
write(mf,*) " Here I set the initial ray of the tracking"
write(mf,*) " "


xcn=(c_n%A_T**(-1)).o.cray


write(mf,*) " "
write(mf,*) " Evaluate the resonance and tune shift"
write(mf,*) " "

bet=c_n0%v(1).o.xcn

dtune=(c_n0%v(2).o.xcn)/twopi

write(mf,*) " Numerical value of beta ",bet

write(mf,*) " Numerical value of beta*beta*",abs(bet)**2

write(mf,*) " Numerical spin tune shift ",dtune 

n=200000
r=0.001d0
dr=2*r/n   ! rate of change of the approach to the resonance


rate=abs(bet)**2/dr/4.d0
write(6,*)"formula with complex DA", 2*exp(-rate)-1
write(mf,*)"formula with complex DA", 2*exp(-rate)-1


call kanalnummer(mft,"spin_barber_rf.dat")
 

close(mf)

!!! These are in initial coordinates
fac=1.7d0
X0=0.D0
X0(1)=0.14D0*fac
X0(3)=0.01D0*fac
X0(5)=1.d0



ag0=(als%start%ag*gamma0-r)/gamma0



p=>als%start
do k=1,als%n
  p%ag=ag0
  p=>p%next
 enddo
XS0=CLOSED
state=default0+only_4d0+SPIN0
    CALL FIND_ORBIT_probe_spin(ALS,XS0,STATE,c_1d_5,FIBRE1=1,theta0=theta0)
    
xs0%ac%x(1)=1.d0
xs0%ac%x(2)=0.d0
    
XS0%X(1)=XS0%X(1)+X0(1)
XS0%X(3)=XS0%X(3)+X0(3)
XS0%s(1)%x(1)=0.d0
XS0%s(1)%x(2)=1.d0
XS0%s(1)%x(3)=0.d0

state=default0+only_4d0+SPIN0+modulation0
 
do i=1,n   !+5000
 
 if(i<=n) then
 ag=(ag0*gamma0+i*dr)/gamma0
   p=>als%start
 do k=1,als%n
  p%ag=ag
  p=>p%next
 enddo
 endif
 
 CALL TRACK_PROBE(ALS,XS0,STATE,FIBRE1=1)

 IF(MOD(I,10000)==0) then
  write(6,*) i*100.d0/n ," % done"
 endif
 IF(MOD(I,100)==1) then
  write(mft,'(1x,i6,2(1x,E15.7))') i,(i-n/2)*dr,xs0%s(1)%x(2)
 endif
  IF(.NOT.CHECK_STABLE) then 
   pause 999
   stop
  endif
enddo

close(mft)



    call KILL(ray8)
    call KILL(id_S,one_turn_S,A_S,A_f,A_spin,A_l,A_nl,NORMALSPIN)
    call KILL(nf_S)
    call KILL(xs)
    call KILL(theta)
    call KILL(N_AXIS_res)




 call ptc_end


end program ptc_geometry

subroutine  build_ALS(ALS,mis)
use madx_ptc_module
use pointer_lattice
implicit none

type(layout), target :: ALS

real(dp) :: alpha,lbend, cut, ksd, ksf 
type(fibre)  L1,L2,L3,L4,L5,L6,L7,L8,L9,L10
type(fibre)  L11,L12,L13,L14,L15,L16,L17,L18,L19,L20
type(fibre)  L21,L22,L23,L24,L25,L26,L27,L27A,L27B,L27C,L27D,DS
 type(fibre)  QF1,QF2,QD1,QD2,QFA1,QFA2,sf,sd,cav,bend,vc5,bend1
type(layout) :: sfline,sdline,sup1,supb
logical(lp) mis
!-----------------------------------

call make_states(.true.)
exact_model = .false.
!default = default + nocavity  
call update_states
madlength = .false.


call set_mad(energy = 1.5d0, method = 6, step = 3)

madkind2 = drift_kick_drift


  L1  = drift("L1 ",  2.832695d0)
  L2  = drift("L2 ",  0.45698d0)
  L3  = drift("L3 ",  0.08902d0)
  L4  = drift("L4 ",  0.2155d0)
  L5  = drift("L5 ",  0.219d0)
  L6  = drift("L6 ",  0.107078d0)
  L7  = drift("L7 ",  0.105716d0)
  L8  = drift("L8 ",  0.135904d0)
  L9  = drift("L9 ",  0.2156993d0)
  L10 = drift("L10",  0.089084d0)
   L11= drift("L11",  0.235416d0)
   L12= drift("L12",  0.1245d0)
   L13= drift("L13",  0.511844d0)
   L14= drift("L14",  0.1788541d0)
   L15= drift("L15",  0.1788483d0)
   L16= drift("L16",  0.511849d0)
   L17= drift("L17",  0.1245d0)
   L18= drift("L18",  0.235405d0)
   L19= drift("L19",  0.089095d0)
   L20= drift("L20",  0.2157007d0)
   L21= drift("L21",  0.177716d0)
   L22= drift("L22",  0.170981d0)
   L23= drift("L23",  0.218997d0)
 L24 = drift ("L24",  0.215503d0)
 L25 = drift ("L25",  0.0890187d0)
 L26 = drift ("L26",  0.45698d0)
 L27 = drift ("L27",  2.832696d0)
 L27a  = drift (" L27a",  0.8596d0)
 L27b  = drift (" L27b",  0.1524d0)
 L27c  = drift (" L27c",  0.04445d0)
 L27d  = drift (" L27d",  1.776246d0)
 ds  = drift (" DS  ", 0.1015d0)

  QF1 = QUADRUPOLE(" QF1 ",0.344D0, K1= 2.2474D0+6.447435260914397D-03)
  QF2 = QUADRUPOLE(" QF2 ",0.344D0, K1= 2.2474D0)
  QD1 = QUADRUPOLE(" QD1 ",0.187D0, K1= -2.3368D0-2.593018157427161D-02); 
  QD2 = QUADRUPOLE(" QD2 ",0.187D0, K1= -2.3368D0);  
  QFA1= QUADRUPOLE(" QFA1",0.448D0, K1= 2.8856D0);  
  QFA2= QUADRUPOLE(" QFA2",0.448D0, K1= 2.8856D0);  

ksf=(-41.67478927130080+0.3392376315938252)*2.0;ksd= (56.36083889436033-0.1043679358857811)*2.d0;
   sf=sextupole ("sf",2.d0*0.1015d0, K2= ksf);
   sd= sextupole("sd", 2.d0*0.1015d0, K2= ksd);

 VC5=marker("vc5");
ALPHA=0.17453292519943295769236907684886;
 
LBEND=0.86621;
 
 
BEND = RBEND("BEND", LBEND, ANGLE=ALPHA).q.(-0.778741d0)
BEND1 = RBEND("BEND1", LBEND, ANGLE=ALPHA).q.(-0.778741d0)
 
CAV=RFCAVITY("CAV",L=0.0000d0,VOLT=-1.0d0,REV_FREQ=500.0d6)


sfline=1*sf;
sdline=1*sd;

SUP1=L1+L2+L3+QF1+VC5+L4+L5+QD1+L6+L7+L8+VC5+BEND+VC5+L9+sfline+L10+&
           L11+QFA1+L12+sdline+L13+ &
           L14+BEND+L15+L16+sdline+L17+ &
           QFA2+L18+L19+sfline+L20+BEND+L21+&
           L22+QD2+L23+L24+QF2+L25+ &
           L26+VC5+L27;

SUPb=L1+L2+L3+QF1+VC5+L4+L5+QD1+L6+L7+L8+VC5+BEND+VC5+L9+sfline+L10+&
           L11+QFA1+L12+sdline+L13+ &
           L14+BEND+L15+L16+sdline+L17+ &
           QFA2+L18+L19+sfline+L20+BEND1+L21+&
           L22+QD2+L23+L24+QF2+L25+ &
           L26+VC5+L27;

ALS = 11*sup1+supb+cav;
 
ALS = .ring.ALS

call survey(ALS)


if(mis) then
 sig=1.d-5
 cut=4.d0
 call MESS_UP_ALIGNMENT(ALS,SIG,cut)
endif
end subroutine build_ALS
